---
title: "Starbility"
author: "Julia Ma"
date: "2024-08-08"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)

library(haven)
library(tidyverse)
library(starbility)
library(lfe)

```

## Load data
setwd("/Users/ag734/Dropbox/JPE_Replication/Data/Final")

```{r load data}
data <- read_dta("/Users/ag734/Dropbox/JPE_Replication/Data/Final/linesectionlevel_starbility.dta")

```

## Starbility

```{r regression}

# Categorical variables
data$dependency <- as.factor(data$dependency)
data$line_section <- as.factor(data$line_section)
data$Shift <- as.factor(data$Shift)
data$doy <- as.factor(data$doy)
data$line_sec_team <- as.factor(data$line_section_team)
data$line_team <- as.factor(data$line_team)


# Interaction terms
data$mixed_dep0 <- with(data, ifelse(dependency == 0, mixed, 0))
data$mixed_dep1 <- with(data, ifelse(dependency == 1, mixed, 0))

data$tenure_skilled_dep0 <- with(data, ifelse(dependency == 0, tenure_skilled_std, 0))
data$tenure_skilled_dep1 <- with(data, ifelse(dependency == 1, tenure_skilled_std, 0))
data$tenure_unskilled_dep0 <- with(data, ifelse(dependency == 0, tenure_unskilled_std, 0))
data$tenure_unskilled_dep1 <- with(data, ifelse(dependency == 1, tenure_unskilled_std, 0))
data$school_years_dep0 <- with(data, ifelse(dependency == 0, school_years_std, 0))
data$school_years_dep1 <- with(data, ifelse(dependency == 1, school_years_std, 0))

data$mixed_count_workers_std <- with(data, mixed * count_workers_std)
data$mixed_hs_hindu_std <- with(data, mixed * share_hskilled_hindus_std)
data$mixed_hs_muslim_std <- with(data, mixed * share_hskilled_muslims_std)

data$mixed_agestd <- with(data, mixed * age_std)
data$mixed_tenures <- with(data, mixed * tenure_skilled_std)
data$mixed_tenureus <- with(data, mixed * tenure_unskilled_std)
data$mixedXcountw_skilled <- with(data,mixed*count_high_skilled_std)
data$mixedXcountw_unskilled <- with(data,mixed*count_unskilled_std)



# Regression with fixed effects and clustering
model <- felm(
  rating ~ mixed_dep0 + mixed_dep1 + mixed_count_workers_std +
  tenure_skilled_dep0 + tenure_skilled_dep1 +
  tenure_unskilled_dep0 + tenure_unskilled_dep1 +
  school_years_dep0 + school_years_dep1 +
  mixed_hs_hindu_std + mixed_hs_muslim_std +  mixed_tenures + mixed_tenureus + age_std + mixed_agestd + line_section + Shift | doy | 0 | line_sec_team,
  data = data
)

# Summary
summary(model)
```

```{r controls}

# Define controls
#perm_controls <- c(
 # "Mixed × Group Size" = "mixed_count_workers_std",
  #"LD × Tenure (Skilled)" = "tenure_skilled_dep0",
  #"HD × Tenure (Skilled)" = "tenure_skilled_dep1",
  #"LD × Tenure (Unskilled)" = "tenure_unskilled_dep0",
  #"HD × Tenure (Unskilled)" = "tenure_unskilled_dep1",
  #"LD × Years of School" = "school_years_dep0",
  #"HD × Years of School" = "school_years_dep1",
  #"Mixed × Share Skilled (Hindu)" = "mixed_hs_hindu_std",
  #"Mixed × Share Skilled (Muslim)" = "mixed_hs_muslim_std"
#)

#Define controls
perm_controls <- c(
  "Tenure" = "tenure_skilled_dep0 + tenure_skilled_dep1 + tenure_unskilled_dep0 + tenure_unskilled_dep1 +  mixed_tenures + mixed_tenureus",
  "Schooling" = "school_years_dep0 + school_years_dep1",
  "Age" = "mixed_agestd",
  "Team size & Share Semi-skilled/Operator (by religion)" = "mixed_hs_hindu_std + mixed_hs_muslim_std + mixedXcountw_skilled + mixedXcountw_unskilled"
)



base_controls <- c(
  "Day" = "doy",
  "Shift" = "Shift",
  "Line × Section" = "line_section"
)

cluster_controls <- "line_sec_team"
```


```{r starbility plot}

# Stability plot for mixed_dep0
plot_dep0 <- stability_plot(
  data = data,
  lhs = "rating",
  rhs = "mixed_dep0",  # Main independent variable for dep = 0
  perm = perm_controls,
  base = base_controls,
  cluster = cluster_controls,
  rel_height = 0.6,
  error_geom = 'ribbon'
)

# Print
print(plot_dep0)

# Stability plot for mixed_dep1
plot_dep1 <- stability_plot(
  data = data,
  lhs = "rating",
  rhs = "mixed_dep1",  # Main independent variable for dep = 1
  perm = perm_controls,
  base = base_controls,
  cluster = cluster_controls,
  rel_height = 0.6,
  error_geom = 'ribbon'
)

# Print
print(plot_dep1)

```